module interface_d

  interface
     subroutine file_read(coef_max, dim)
!  subroutine file_read(coef_max, b_vec)
    implicit none
    real,intent(inout), allocatable, dimension(:,:) :: coef_max
    !    real, intent(inout), allocatable, dimension(:) :: b_vec
    integer, intent(inout) :: dim
  end subroutine file_read
end interface
end module interface_d

!--------------------------------------------------------------------------
program test
  use interface_d
  implicit none
  real, allocatable, dimension(:,:) :: coef_max
  !  real, allocatable,dimension(:) :: b_vec
  integer :: dim
  
  call file_read(coef_max, dim)
  call ga_eli(coef_max, dim)
  
  deallocate(coef_max)
!  deallocate(b_vec)

end program test
!---------------------------------------------------------------------------     
subroutine file_read(coef_max, dim)
!subroutine file_read(coef_max, b_vec)
  implicit none
  real,intent(inout), allocatable, dimension(:,:) :: coef_max
  !  real, intent(inout), allocatable, dimension(:) :: b_vec
  integer, intent(inout) :: dim
  integer ::  k, status
  character(len=20) :: fname

  write(*,*) 'input file:'
  read(*,*) fname
  open(unit=7, file=fname, status='old', action='read', iostat=status)

  openif: if(status /= 0) then
     write(*,*) 'open file error'
  else
     read(7,*) dim
     allocate(coef_max(dim,dim+1))
 !    allocate(b_vec(dim))
     
     read_coef_max: do k = 1,dim
        read(7,*) coef_max(k,:)
     end do read_coef_max

  !   read(7,*) b_vec
  end if openif
  call pnt(coef_max, dim, dim+1)
  
end subroutine file_read
!------------------------------------------------------------------------------------
subroutine ga_eli(coef_max, dim)
  implicit none
  integer, intent(in) :: dim
  real, intent(inout), dimension(dim,dim+1) :: coef_max
  real, dimension(dim) :: solution
  integer :: k, k2, k3, index
  real :: tm, t, b, ma

  do k = 1, dim-1
     
     ma = abs(coef_max(k,k))
     find_max: do k2 = k+1,dim
        if(abs(coef_max(k2, k)) > ma) then
           ma = abs(coef_max(k2, k));
           index = k2
        end if
     end do find_max
!     write(*,*) ' find_max'
!    write(*,*) ma, index
     
     change_line: do k2 = k, dim+1
        t = coef_max(k, k2)
        coef_max(k, k2) = coef_max(index, k2)
        coef_max(index, k2) = t
     end do change_line
!     call pnt(coef_max, dim, dim+1)
     
     eliminate: do k2 = k+1, dim
        b = coef_max(k2, k) / coef_max(k, k)
        write(*,*) b
        sub: do k3 = k+1, dim+1
           coef_max(k2, k3) = coef_max(k2, k3) - b*coef_max(k,k3)
        end do sub
     end do eliminate
!    call pnt(coef_max, dim, dim+1)
  end do

  solution(dim) = coef_max(dim,dim+1) / coef_max(dim,dim)
  back: do k = dim-1,1,-1
     tm = coef_max(k,dim+1)
     do k2 = dim,k+1, -1
        tm = tm - solution(k2) * coef_max(k,k2)
     end do
     solution(k) = tm / coef_max(k,K)
  end do back

  write(*,*) solution
        

end subroutine ga_eli


subroutine pnt(a, row, col)
  integer, intent(in) :: row, col
  real, intent(in), dimension(row,col) :: a
  integer :: k
  write(*,*) 
  do k = 1, row
     write(*,*) a(k,:)
  end do
end subroutine pnt
